1) Problem 14.3

#import data
data <- read_excel("Chapter14_exercises_data.xls", sheet = "Exercise 3, 4, 10")
sp500_ts <- ts(data$`Adj Close`, start = 2000, freq = 252)
plot(sp500_ts)

There is noticeable high volatility around the early 2000s and 2008, corresponding to the dotcom bubble and the financial crisis, respectively. These periods were marked by significant instability and substantial fluctuations in stock prices. In contrast, there were also calmer periods of lower volatility, such as the stretch from 2003 to 2007.

#get daily returns
returns <- diff(log(sp500_ts))
tsdisplay(returns, lag.max = 40)

returns_sq <- returns^2
tsdisplay(returns_sq, lag.max = 40)

The returns data seem to exhibit white noise, with the acf/pacf plots not really showing decay or significant autocorrelation. The amplitude of 0.05 is quite low as well. After squaring the data and looking at the acf/pacf plots, there is decay and high persistence in the ACF, with spikes cutting off at lag 14 in the PACF. I will try fitting an ARCH(14).

model=ugarchspec(
  variance.model = list(model = "sGARCH", garchOrder = c(14, 0)),
  mean.model = list(armaOrder = c(2, 2), include.mean = TRUE),
  distribution.model = "sstd")

#sanity check: explore the model parameters
model@model$pars
##          Level Fixed Include Estimate LB UB
## mu           0     0       1        1 NA NA
## ar1          0     0       1        1 NA NA
## ar2          0     0       1        1 NA NA
## ma1          0     0       1        1 NA NA
## ma2          0     0       1        1 NA NA
## arfima       0     0       0        0 NA NA
## archm        0     0       0        0 NA NA
## mxreg        0     0       0        0 NA NA
## omega        0     0       1        1 NA NA
## alpha1       0     0       1        1 NA NA
## alpha2       0     0       1        1 NA NA
## alpha3       0     0       1        1 NA NA
## alpha4       0     0       1        1 NA NA
## alpha5       0     0       1        1 NA NA
## alpha6       0     0       1        1 NA NA
## alpha7       0     0       1        1 NA NA
## alpha8       0     0       1        1 NA NA
## alpha9       0     0       1        1 NA NA
## alpha10      0     0       1        1 NA NA
## alpha11      0     0       1        1 NA NA
## alpha12      0     0       1        1 NA NA
## alpha13      0     0       1        1 NA NA
## alpha14      0     0       1        1 NA NA
## beta         0     0       0        0 NA NA
## gamma        0     0       0        0 NA NA
## eta1         0     0       0        0 NA NA
## eta2         0     0       0        0 NA NA
## delta        0     0       0        0 NA NA
## lambda       0     0       0        0 NA NA
## vxreg        0     0       0        0 NA NA
## skew         0     0       1        1 NA NA
## shape        0     0       1        1 NA NA
## ghlambda     0     0       0        0 NA NA
## xi           0     0       0        0 NA NA
#fit the model to the data
modelfit=ugarchfit(spec=model,data=returns)
modelfit
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(14,0)
## Mean Model   : ARFIMA(2,0,2)
## Distribution : sstd 
## 
## Optimal Parameters
## ------------------------------------
##          Estimate  Std. Error   t value Pr(>|t|)
## mu       0.000424    0.000115  3.694424 0.000220
## ar1      0.028714    0.333601  0.086074 0.931407
## ar2      0.379502    0.217792  1.742496 0.081422
## ma1     -0.106715    0.330267 -0.323117 0.746606
## ma2     -0.446759    0.238626 -1.872216 0.061177
## omega    0.000017    0.000003  6.637722 0.000000
## alpha1   0.000000    0.014376  0.000000 1.000000
## alpha2   0.124089    0.025264  4.911631 0.000001
## alpha3   0.090595    0.023149  3.913513 0.000091
## alpha4   0.092997    0.024990  3.721423 0.000198
## alpha5   0.104300    0.024761  4.212223 0.000025
## alpha6   0.061165    0.022886  2.672616 0.007526
## alpha7   0.068429    0.023645  2.894055 0.003803
## alpha8   0.096837    0.025670  3.772333 0.000162
## alpha9   0.052590    0.021073  2.495599 0.012574
## alpha10  0.072525    0.024225  2.993823 0.002755
## alpha11  0.072022    0.022884  3.147225 0.001648
## alpha12  0.022210    0.020987  1.058290 0.289923
## alpha13  0.036099    0.019535  1.847881 0.064620
## alpha14  0.023143    0.019559  1.183250 0.236710
## skew     0.870222    0.021124 41.195553 0.000000
## shape    8.630207    1.297672  6.650531 0.000000
## 
## Robust Standard Errors:
##          Estimate  Std. Error  t value Pr(>|t|)
## mu       0.000424    0.000129  3.29899 0.000970
## ar1      0.028714    0.216027  0.13292 0.894256
## ar2      0.379502    0.114688  3.30898 0.000936
## ma1     -0.106715    0.213577 -0.49966 0.617316
## ma2     -0.446759    0.126705 -3.52599 0.000422
## omega    0.000017    0.000002  7.76844 0.000000
## alpha1   0.000000    0.017545  0.00000 1.000000
## alpha2   0.124089    0.029129  4.25998 0.000020
## alpha3   0.090595    0.021727  4.16974 0.000030
## alpha4   0.092997    0.025845  3.59828 0.000320
## alpha5   0.104300    0.023882  4.36731 0.000013
## alpha6   0.061165    0.023226  2.63348 0.008451
## alpha7   0.068429    0.024684  2.77223 0.005567
## alpha8   0.096837    0.027525  3.51820 0.000434
## alpha9   0.052590    0.020928  2.51295 0.011973
## alpha10  0.072525    0.025149  2.88376 0.003930
## alpha11  0.072022    0.022984  3.13356 0.001727
## alpha12  0.022210    0.023142  0.95973 0.337192
## alpha13  0.036099    0.018060  1.99882 0.045628
## alpha14  0.023143    0.021759  1.06361 0.287507
## skew     0.870222    0.020989 41.46165 0.000000
## shape    8.630207    1.277840  6.75375 0.000000
## 
## LogLikelihood : 10562.67 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.2854
## Bayes        -6.2453
## Shibata      -6.2855
## Hannan-Quinn -6.2711
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic   p-value
## Lag[1]                       1.738 1.874e-01
## Lag[2*(p+q)+(p+q)-1][11]    11.264 1.024e-12
## Lag[4*(p+q)+(p+q)-1][19]    17.604 3.533e-03
## d.o.f=4
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                      0.2877  0.5917
## Lag[2*(p+q)+(p+q)-1][41]   13.3219  0.9431
## Lag[4*(p+q)+(p+q)-1][69]   23.4313  0.9680
## d.o.f=14
## 
## Weighted ARCH LM Tests
## ------------------------------------
##              Statistic Shape Scale P-Value
## ARCH Lag[15]     1.900 0.500 2.000  0.1680
## ARCH Lag[17]     2.027 1.496 1.887  0.5408
## ARCH Lag[19]     2.843 2.483 1.802  0.6715
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  no.parameters>20 (not available)
## Individual Statistics:               
## mu      0.71695
## ar1     0.08409
## ar2     0.04561
## ma1     0.09545
## ma2     0.04949
## omega   0.39837
## alpha1  0.38957
## alpha2  0.32477
## alpha3  0.22171
## alpha4  0.06361
## alpha5  0.30666
## alpha6  0.06821
## alpha7  0.19388
## alpha8  0.42707
## alpha9  0.06410
## alpha10 0.06916
## alpha11 0.23119
## alpha12 0.09882
## alpha13 0.13785
## alpha14 0.09918
## skew    0.31054
## shape   0.47752
## 
## Asymptotic Critical Values (10% 5% 1%)
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value      prob sig
## Sign Bias           1.8962 5.803e-02   *
## Negative Sign Bias  0.6498 5.159e-01    
## Positive Sign Bias  2.7377 6.220e-03 ***
## Joint Effect       30.4067 1.133e-06 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     49.09    0.0001782
## 2    30     56.61    0.0016004
## 3    40     65.99    0.0044346
## 4    50     69.67    0.0276760
## 
## 
## Elapsed time : 6.710096
plot(modelfit, which = 10)

plot(modelfit, which = 11)

Looking at the ACF of squared standardized residuals, it seems that the values stay within the bands, meaning dynamics have been wiped out. This model is good for fitting the data.

#After comparing AIC/BICs for various GARCH combinations, I decided to go with an GARCH(2,1) model, which had the lowest value.
model1=ugarchspec(
  variance.model = list(model = "sGARCH", garchOrder = c(2, 1)),
  mean.model = list(armaOrder = c(2, 2), include.mean = TRUE),
  distribution.model = "sstd")

#sanity check: explore the model parameters
model1@model$pars
##          Level Fixed Include Estimate LB UB
## mu           0     0       1        1 NA NA
## ar1          0     0       1        1 NA NA
## ar2          0     0       1        1 NA NA
## ma1          0     0       1        1 NA NA
## ma2          0     0       1        1 NA NA
## arfima       0     0       0        0 NA NA
## archm        0     0       0        0 NA NA
## mxreg        0     0       0        0 NA NA
## omega        0     0       1        1 NA NA
## alpha1       0     0       1        1 NA NA
## alpha2       0     0       1        1 NA NA
## beta1        0     0       1        1 NA NA
## gamma        0     0       0        0 NA NA
## eta1         0     0       0        0 NA NA
## eta2         0     0       0        0 NA NA
## delta        0     0       0        0 NA NA
## lambda       0     0       0        0 NA NA
## vxreg        0     0       0        0 NA NA
## skew         0     0       1        1 NA NA
## shape        0     0       1        1 NA NA
## ghlambda     0     0       0        0 NA NA
## xi           0     0       0        0 NA NA
#fit the model to the data
modelfit1=ugarchfit(spec=model1,data=returns)
modelfit1
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(2,1)
## Mean Model   : ARFIMA(2,0,2)
## Distribution : sstd 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.000414    0.000101  4.082633 0.000045
## ar1     0.075122    0.347045  0.216460 0.828629
## ar2     0.364454    0.232917  1.564738 0.117644
## ma1    -0.155823    0.344871 -0.451829 0.651392
## ma2    -0.426434    0.257215 -1.657887 0.097340
## omega   0.000002    0.000006  0.293893 0.768839
## alpha1  0.000000    0.014750  0.000032 0.999974
## alpha2  0.114612    0.089699  1.277738 0.201342
## beta1   0.876121    0.087667  9.993787 0.000000
## skew    0.866983    0.012883 67.299125 0.000000
## shape   8.626200    2.659506  3.243535 0.001181
## 
## Robust Standard Errors:
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.000414    0.000879  0.470837 0.637757
## ar1     0.075122    0.156581  0.479760 0.631398
## ar2     0.364454    0.902869  0.403662 0.686461
## ma1    -0.155823    0.154911 -1.005884 0.314471
## ma2    -0.426434    0.902482 -0.472512 0.636561
## omega   0.000002    0.000094  0.017788 0.985808
## alpha1  0.000000    0.075751  0.000006 0.999995
## alpha2  0.114612    1.516198  0.075592 0.939744
## beta1   0.876121    1.457721  0.601021 0.547826
## skew    0.866983    0.276464  3.135972 0.001713
## shape   8.626200   48.751382  0.176943 0.859553
## 
## LogLikelihood : 10563.23 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -6.2923
## Bayes        -6.2723
## Shibata      -6.2923
## Hannan-Quinn -6.2851
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic   p-value
## Lag[1]                       1.598 2.061e-01
## Lag[2*(p+q)+(p+q)-1][11]    11.520 1.323e-13
## Lag[4*(p+q)+(p+q)-1][19]    18.089 2.305e-03
## d.o.f=4
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                    0.005113  0.9430
## Lag[2*(p+q)+(p+q)-1][8]   2.196073  0.8277
## Lag[4*(p+q)+(p+q)-1][14]  3.905436  0.8891
## d.o.f=3
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[4]  0.002388 0.500 2.000  0.9610
## ARCH Lag[6]  1.155457 1.461 1.711  0.7036
## ARCH Lag[8]  1.802242 2.368 1.583  0.7812
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  97.0513
## Individual Statistics:               
## mu      0.73947
## ar1     0.09432
## ar2     0.05534
## ma1     0.10374
## ma2     0.06097
## omega  20.67564
## alpha1  0.20342
## alpha2  0.18104
## beta1   0.17380
## skew    0.38641
## shape   0.49342
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          2.49 2.75 3.27
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value      prob sig
## Sign Bias           1.9387 5.262e-02   *
## Negative Sign Bias  0.5138 6.074e-01    
## Positive Sign Bias  2.7341 6.288e-03 ***
## Joint Effect       29.9240 1.432e-06 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     41.79     0.001892
## 2    30     54.37     0.002931
## 3    40     56.81     0.032544
## 4    50     74.62     0.010622
## 
## 
## Elapsed time : 1.991244
plot(modelfit1, which = 10)

plot(modelfit1, which = 11)

The GARCH model fits equally well since there are generally no spikes. It also has a slightly lower AIC/BIC. I will prefer the GARCH model due to having less parameters compared to the ARCH model.

2) Problem 14.4

#forecast
modelfor = ugarchforecast(modelfit1, data = NULL, n.ahead = 2, n.roll = 0, out.sample = 0)
plot(modelfor, which = 1)

#get forecasted volatility
forecasted_volatility <- sigma(modelfor)

#get 95% prediction intervals for 1 and 2 step ahead forecasts
lower_1 <- -1.96 * forecasted_volatility[1]
upper_1 <- 1.96 * forecasted_volatility[1]
lower_2 <- -1.96 * forecasted_volatility[2]
upper_2 <- 1.96 * forecasted_volatility[2]

#print 95% prediction intervals
cat("One-Step Ahead - Lower Bound:", lower_1, "Upper Bound:", upper_1, "\n")
## One-Step Ahead - Lower Bound: -0.01713364 Upper Bound: 0.01713364
cat("Two-Step Ahead - Lower Bound:", lower_2, "Upper Bound:", upper_2, "\n")
## Two-Step Ahead - Lower Bound: -0.01634723 Upper Bound: 0.01634723

3) Problem 14.5

cpi <- read_excel("Chapter14_exercises_data.xls", sheet = "Exercise 5a")
gdp <- read_excel("Chapter14_exercises_data.xls", sheet = "Exercise 5b")
inf <- cpi$`INFLATION RATE`[2:nrow(cpi)]
gdp_grow <- gdp$`GDP GROWTH`[2:nrow(gdp)]
mean(inf)
## [1] 0.3009486
mean(gdp_grow)
## [1] 0.7842689
plot(inf, type="l")

plot(gdp_grow, type = "l")

adf.test(cpi$`INFLATION RATE`[2:nrow(cpi)])
## 
##  Augmented Dickey-Fuller Test
## 
## data:  cpi$`INFLATION RATE`[2:nrow(cpi)]
## Dickey-Fuller = -3.9489, Lag order = 9, p-value = 0.01155
## alternative hypothesis: stationary
adf.test(gdp$`GDP GROWTH`[2:nrow(gdp)])
## Warning in adf.test(gdp$`GDP GROWTH`[2:nrow(gdp)]): p-value smaller than
## printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  gdp$`GDP GROWTH`[2:nrow(gdp)]
## Dickey-Fuller = -6.4572, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary
auto.arima(cpi$"INFLATION RATE"[2:nrow(cpi)]) #ARIMA order (4,1,1)
## Series: cpi$"INFLATION RATE"[2:nrow(cpi)] 
## ARIMA(4,1,1) 
## 
## Coefficients:
##          ar1      ar2     ar3      ar4      ma1
##       0.2872  -0.0455  0.0071  -0.0227  -0.8607
## s.e.  0.0715   0.0509  0.0495   0.0525   0.0626
## 
## sigma^2 = 0.078:  log likelihood = -111.72
## AIC=235.45   AICc=235.55   BIC=263.5
auto.arima(gdp$"GDP GROWTH"[2:nrow(gdp)]) #ARIMA order (2,0,2)
## Series: gdp$"GDP GROWTH"[2:nrow(gdp)] 
## ARIMA(2,0,2) with non-zero mean 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2    mean
##       1.2589  -0.6531  -0.9388  0.4793  0.7814
## s.e.  0.1742   0.1593   0.2075  0.1654  0.0763
## 
## sigma^2 = 0.8324:  log likelihood = -348
## AIC=708   AICc=708.33   BIC=729.46
diff_inf<- diff(inf)
#Model for inflation rate
model <- ugarchspec(
  variance.model = list(model = "sGARCH", garchOrder = c(3, 3)),
  mean.model = list(armaOrder = c(4, 1), include.mean = TRUE),
  distribution.model = "sstd"
)

# Sanity check: explore the model parameters
print(model@model$pars)
##          Level Fixed Include Estimate LB UB
## mu           0     0       1        1 NA NA
## ar1          0     0       1        1 NA NA
## ar2          0     0       1        1 NA NA
## ar3          0     0       1        1 NA NA
## ar4          0     0       1        1 NA NA
## ma1          0     0       1        1 NA NA
## arfima       0     0       0        0 NA NA
## archm        0     0       0        0 NA NA
## mxreg        0     0       0        0 NA NA
## omega        0     0       1        1 NA NA
## alpha1       0     0       1        1 NA NA
## alpha2       0     0       1        1 NA NA
## alpha3       0     0       1        1 NA NA
## beta1        0     0       1        1 NA NA
## beta2        0     0       1        1 NA NA
## beta3        0     0       1        1 NA NA
## gamma        0     0       0        0 NA NA
## eta1         0     0       0        0 NA NA
## eta2         0     0       0        0 NA NA
## delta        0     0       0        0 NA NA
## lambda       0     0       0        0 NA NA
## vxreg        0     0       0        0 NA NA
## skew         0     0       1        1 NA NA
## shape        0     0       1        1 NA NA
## ghlambda     0     0       0        0 NA NA
## xi           0     0       0        0 NA NA
# Fit the model to the data
modelfit1 <- ugarchfit(spec = model, data = diff_inf)
print(modelfit1)
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(3,3)
## Mean Model   : ARFIMA(4,0,1)
## Distribution : sstd 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error    t value Pr(>|t|)
## mu      0.001238    0.001359   0.911414 0.362077
## ar1     0.196110    0.060187   3.258320 0.001121
## ar2    -0.004106    0.049146  -0.083547 0.933417
## ar3    -0.042455    0.046840  -0.906382 0.364734
## ar4    -0.019913    0.043901  -0.453588 0.650125
## ma1    -0.838828    0.043947 -19.087256 0.000000
## omega   0.008984    0.011889   0.755648 0.449860
## alpha1  0.197500    0.072165   2.736772 0.006205
## alpha2  0.168374    0.344144   0.489255 0.624661
## alpha3  0.124121    0.233845   0.530782 0.595570
## beta1   0.000000    1.771586   0.000000 1.000000
## beta2   0.320327    0.786156   0.407460 0.683670
## beta3   0.120003    0.262842   0.456561 0.647987
## skew    1.086276    0.049806  21.809933 0.000000
## shape   4.698211    0.676905   6.940723 0.000000
## 
## Robust Standard Errors:
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.001238    0.001565  0.791306 0.428766
## ar1     0.196110    0.102226  1.918398 0.055061
## ar2    -0.004106    0.077936 -0.052684 0.957984
## ar3    -0.042455    0.076323 -0.556255 0.578036
## ar4    -0.019913    0.082114 -0.242504 0.808390
## ma1    -0.838828    0.116500 -7.200259 0.000000
## omega   0.008984    0.051475  0.174529 0.861450
## alpha1  0.197500    0.103260  1.912642 0.055794
## alpha2  0.168374    1.496250  0.112531 0.910403
## alpha3  0.124121    1.053896  0.117773 0.906247
## beta1   0.000000    7.290230  0.000000 1.000000
## beta2   0.320327    3.509238  0.091281 0.927269
## beta3   0.120003    0.884758  0.135634 0.892111
## skew    1.086276    0.084951 12.787098 0.000000
## shape   4.698211    1.629176  2.883795 0.003929
## 
## LogLikelihood : 51.65621 
## 
## Information Criteria
## ------------------------------------
##                        
## Akaike       -0.0924495
## Bayes        -0.0040039
## Shibata      -0.0931475
## Hannan-Quinn -0.0584586
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic   p-value
## Lag[1]                   6.373e-05 9.936e-01
## Lag[2*(p+q)+(p+q)-1][14] 1.195e+01 9.639e-11
## Lag[4*(p+q)+(p+q)-1][24] 2.828e+01 3.063e-06
## d.o.f=5
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                          statistic   p-value
## Lag[1]                       11.36 0.0007519
## Lag[2*(p+q)+(p+q)-1][17]     17.38 0.0219348
## Lag[4*(p+q)+(p+q)-1][29]     19.65 0.1496115
## d.o.f=6
## 
## Weighted ARCH LM Tests
## ------------------------------------
##              Statistic Shape Scale P-Value
## ARCH Lag[7]    0.08025 0.500 2.000  0.7770
## ARCH Lag[9]    0.84418 1.485 1.796  0.8114
## ARCH Lag[11]   0.90917 2.440 1.677  0.9505
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  4.932
## Individual Statistics:              
## mu     0.04198
## ar1    0.24350
## ar2    0.98424
## ar3    0.21016
## ar4    0.10888
## ma1    0.35696
## omega  0.13734
## alpha1 0.19043
## alpha2 0.16278
## alpha3 0.33587
## beta1  0.21233
## beta2  0.30283
## beta3  0.28592
## skew   0.06654
## shape  0.26498
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          3.26 3.54 4.07
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value     prob sig
## Sign Bias           1.4630 0.143874    
## Negative Sign Bias  0.4078 0.683530    
## Positive Sign Bias  3.0559 0.002319 ***
## Joint Effect        9.5390 0.022921  **
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     19.79       0.4075
## 2    30     27.74       0.5316
## 3    40     37.62       0.5329
## 4    50     61.92       0.1018
## 
## 
## Elapsed time : 1.085289
plot(modelfit1, which = 10)

plot(modelfit1, which = 11)

modelfor = ugarchforecast(modelfit1, data = NULL, n.ahead = 2, n.roll = 0, out.sample = 0)
plot(modelfor, which = 1)

#get forecasted volatility
forecasted_volatility <- sigma(modelfor)

#One step ahead volatility forecast
lower_1 <- -1.96 * forecasted_volatility[1]
upper_1 <- 1.96 * forecasted_volatility[1]

#print 95% prediction intervals
cat("One-Step Ahead - Lower Bound:", lower_1, "Upper Bound:", upper_1, "\n")
## One-Step Ahead - Lower Bound: -0.7550893 Upper Bound: 0.7550893
#Model for GDP Growth
model_2 <- ugarchspec(
  variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(2, 2), include.mean = TRUE),
  distribution.model = "sstd"
)

# Sanity check: explore the model parameters
print(model_2@model$pars)
##          Level Fixed Include Estimate LB UB
## mu           0     0       1        1 NA NA
## ar1          0     0       1        1 NA NA
## ar2          0     0       1        1 NA NA
## ma1          0     0       1        1 NA NA
## ma2          0     0       1        1 NA NA
## arfima       0     0       0        0 NA NA
## archm        0     0       0        0 NA NA
## mxreg        0     0       0        0 NA NA
## omega        0     0       1        1 NA NA
## alpha1       0     0       1        1 NA NA
## beta1        0     0       1        1 NA NA
## gamma        0     0       0        0 NA NA
## eta1         0     0       0        0 NA NA
## eta2         0     0       0        0 NA NA
## delta        0     0       0        0 NA NA
## lambda       0     0       0        0 NA NA
## vxreg        0     0       0        0 NA NA
## skew         0     0       1        1 NA NA
## shape        0     0       1        1 NA NA
## ghlambda     0     0       0        0 NA NA
## xi           0     0       0        0 NA NA
modelfit_2 <- ugarchfit(spec = model, data = gdp_grow)
print(modelfit_2)
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(3,3)
## Mean Model   : ARFIMA(4,0,1)
## Distribution : sstd 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.768010    0.173963  4.41478 0.000010
## ar1     1.192082    0.229627  5.19139 0.000000
## ar2    -0.071943    0.160867 -0.44722 0.654714
## ar3    -0.342749    0.108333 -3.16386 0.001557
## ar4     0.134246    0.068059  1.97248 0.048554
## ma1    -0.832107    0.246534 -3.37523 0.000738
## omega   0.082043    0.215821  0.38014 0.703839
## alpha1  0.271551    0.118106  2.29921 0.021493
## alpha2  0.212252    0.362164  0.58607 0.557831
## alpha3  0.000000    0.672449  0.00000 1.000000
## beta1   0.000000    0.280433  0.00000 1.000000
## beta2   0.000000    0.657710  0.00000 1.000000
## beta3   0.451791    0.116947  3.86322 0.000112
## skew    1.018403    0.128628  7.91742 0.000000
## shape   7.818062    4.907829  1.59298 0.111165
## 
## Robust Standard Errors:
##         Estimate  Std. Error   t value Pr(>|t|)
## mu      0.768010     1.02994  0.745683  0.45586
## ar1     1.192082     0.14201  8.394126  0.00000
## ar2    -0.071943     0.69750 -0.103145  0.91785
## ar3    -0.342749     0.23745 -1.443454  0.14889
## ar4     0.134246     0.18383  0.730264  0.46523
## ma1    -0.832107     0.68262 -1.218988  0.22285
## omega   0.082043     1.35318  0.060629  0.95165
## alpha1  0.271551     0.26193  1.036711  0.29987
## alpha2  0.212252     2.32039  0.091473  0.92712
## alpha3  0.000000     4.35927  0.000000  1.00000
## beta1   0.000000     1.90540  0.000000  1.00000
## beta2   0.000000     4.61083  0.000000  1.00000
## beta3   0.451791     0.28323  1.595127  0.11068
## skew    1.018403     0.60007  1.697139  0.08967
## shape   7.818062    20.95117  0.373156  0.70903
## 
## LogLikelihood : -319.1496 
## 
## Information Criteria
## ------------------------------------
##                    
## Akaike       2.5314
## Bayes        2.7346
## Shibata      2.5254
## Hannan-Quinn 2.6131
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                      0.5389  0.4629
## Lag[2*(p+q)+(p+q)-1][14]    5.6791  0.9995
## Lag[4*(p+q)+(p+q)-1][24]   11.9004  0.5632
## d.o.f=5
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                      0.0158  0.9000
## Lag[2*(p+q)+(p+q)-1][17]    6.1031  0.7881
## Lag[4*(p+q)+(p+q)-1][29]   11.1602  0.7980
## d.o.f=6
## 
## Weighted ARCH LM Tests
## ------------------------------------
##              Statistic Shape Scale P-Value
## ARCH Lag[7]      1.202 0.500 2.000  0.2730
## ARCH Lag[9]      2.467 1.485 1.796  0.4268
## ARCH Lag[11]     3.157 2.440 1.677  0.5671
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  3.0389
## Individual Statistics:              
## mu     0.19668
## ar1    0.25910
## ar2    0.40438
## ar3    0.42239
## ar4    0.43056
## ma1    0.15904
## omega  0.52367
## alpha1 0.29816
## alpha2 0.38881
## alpha3 0.25425
## beta1  0.61003
## beta2  0.50725
## beta3  0.69710
## skew   0.03317
## shape  0.26153
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          3.26 3.54 4.07
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias          0.90766 0.3649    
## Negative Sign Bias 0.03468 0.9724    
## Positive Sign Bias 0.59490 0.5524    
## Joint Effect       3.39834 0.3342    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     11.00       0.9238
## 2    30     21.45       0.8421
## 3    40     44.18       0.2619
## 4    50     37.89       0.8752
## 
## 
## Elapsed time : 0.4812002
plot(modelfit_2, which = 10)

plot(modelfit_2, which = 11)

modelfor = ugarchforecast(modelfit1, data = NULL, n.ahead = 2, n.roll = 0, out.sample = 0)
plot(modelfor, which = 1)

forecasted_volatility <- sigma(modelfor)

#One step ahead volatility forecast
lower_1 <- -1.96 * forecasted_volatility[1]
upper_1 <- 1.96 * forecasted_volatility[1]

#95% prediction intervals
cat("One-Step Ahead - Lower Bound:", lower_1, "Upper Bound:", upper_1, "\n")
## One-Step Ahead - Lower Bound: -0.7550893 Upper Bound: 0.7550893

4. Problem 12.2

(a)

gasoline_ts <- ts(us_gasoline$Barrels, frequency=52, start = c(1991, 6))
# lambda of Box-Cox transformation
lambda_gasoline <- BoxCox.lambda(gasoline_ts)

# select the number of Fourier pairs
min.AIC <- Inf
K_min.Aic <- 0

for(num in c(1:6)){
  gasoline_tslm <- tslm(
    gasoline_ts ~ trend + fourier(gasoline_ts, K = num),
    lambda = lambda_gasoline
    )
  
  AIC <- CV(gasoline_tslm)["AIC"]
  
  if(AIC < min.AIC){
    min.AIC <- AIC
    K_min.Aic <- num
  }
}

# Make harmonic regression model
gasoline_tslm <- tslm(
  gasoline_ts ~ trend + fourier(gasoline_ts, K = K_min.Aic),
  lambda = lambda_gasoline
  )

autoplot(gasoline_ts) + autolayer(gasoline_tslm$fitted.values)

The dynamic harmonic regression model fits the data better than the harmonic regression model because it is able to capture changes in the seasonality and trend over time.

(b)

# checkresiduals(gasoline_autoarima)
checkresiduals(gasoline_tslm)

## 
##  Breusch-Godfrey test for serial correlation of order up to 104
## 
## data:  Residuals from Linear regression model
## LM test = 932.9, df = 104, p-value < 2.2e-16

The residuals from the dynamic harmonic regression are almost all positive and there are several significant lags in the ACF plot of the residuals, but the autocorrelation is generally quite small. The residuals do not have constant volatility over time. For the harmonic regression model, the residuals are clearly not white noise as they display both trend and changing volatility. The ACFs decay extremely slowly starting from around 0.7.

(c)

We could also try models like ARIMA and ETS, although they might not perform as well on high frequency like weekly data with complex seasonal patterns. We generally used them for monthly or quarterly data.

# Fit ARIMA model
fit_arima <- auto.arima(gasoline_ts)
summary(fit_arima)
## Series: gasoline_ts 
## ARIMA(0,1,1)(1,0,0)[52] 
## 
## Coefficients:
##           ma1    sar1
##       -0.7440  0.2618
## s.e.   0.0206  0.0294
## 
## sigma^2 = 0.07239:  log likelihood = -144.86
## AIC=295.72   AICc=295.74   BIC=311.35
## 
## Training set error measures:
##                       ME     RMSE       MAE         MPE     MAPE      MASE
## Training set 0.003688876 0.268748 0.2067944 -0.04355785 2.480655 0.7176306
##                     ACF1
## Training set -0.01839917
checkresiduals(fit_arima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)(1,0,0)[52]
## Q* = 330.24, df = 102, p-value < 2.2e-16
## 
## Model df: 2.   Total lags used: 104
# Fit ETS model
fit_ets <- ets(gasoline_ts)
## Warning in ets(gasoline_ts): I can't handle data with frequency greater than
## 24. Seasonality will be ignored. Try stlf() if you need seasonal forecasts.
summary(fit_ets)
## ETS(A,N,N) 
## 
## Call:
##  ets(y = gasoline_ts) 
## 
##   Smoothing parameters:
##     alpha = 0.3251 
## 
##   Initial states:
##     l = 6.7151 
## 
##   sigma:  0.2772
## 
##      AIC     AICc      BIC 
## 6299.071 6299.089 6314.705 
## 
## Training set error measures:
##                       ME      RMSE       MAE         MPE     MAPE      MASE
## Training set 0.003778683 0.2770348 0.2132618 -0.04538429 2.560032 0.7400743
##                      ACF1
## Training set -0.002521048
checkresiduals(fit_ets)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,N,N)
## Q* = 455.07, df = 104, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 104
autoplot(gasoline_ts, ylab = "Gas Supply (Weekly)",main= "Other Models") +  autolayer(fitted(fit_arima)) + autolayer(fitted(fit_ets))

5) 12.3

# experiment with retail data.
# create retail data
retail <- read.xlsx("retail.xlsx",
                    sheetIndex = 1,
                    startRow = 2)
retail.ts <- ts(retail[,"A3349873A"], 
                frequency=12, 
                start=c(1982,4))
# fit model
retail_nnetar <- nnetar(retail.ts)
checkresiduals(retail_nnetar)

## 
##  Ljung-Box test
## 
## data:  Residuals from NNAR(3,1,2)[12]
## Q* = 311.51, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24
autoplot(retail.ts) + autolayer(retail_nnetar$fitted)

# Produce forecasts
fc_retail_nnetar <- forecast(retail_nnetar, h = 36)
autoplot(fc_retail_nnetar)

# test accuracy using future data.
retail.new <- read.xlsx("8501011.xlsx", 
                        sheetName = "Data1", 
                        startRow = 10)
retail.new.ts <- ts(retail.new[, "A3349873A"],
                    start = c(1982, 4),
                    frequency = 12)
retail.new.test <- subset(
  retail.new.ts,
  start = length(retail.ts) + 1
  )

# test accuracy using future data.
accuracy(fc_retail_nnetar, retail.new.test)
##                        ME     RMSE      MAE         MPE      MAPE      MASE
## Training set   0.01200432 19.38566 13.77693  -0.7347354  6.078575 0.7276128
## Test set     -51.12171115 69.70154 64.71798 -12.0187699 14.337653 3.4180073
##                   ACF1 Theil's U
## Training set 0.5539058        NA
## Test set     0.6514395  1.434885

Even the the residuals aren’t like white nice and were slightly right skewed, it looks like neural network model fitted well to the data.

# experiment with ibmclose data.
ibmclose_nnetar <- nnetar(ibmclose)
# Produce forecasts
autoplot(ibmclose) + autolayer(ibmclose_nnetar$fitted)

checkresiduals(ibmclose_nnetar)

## 
##  Ljung-Box test
## 
## data:  Residuals from NNAR(1,1)
## Q* = 15.703, df = 10, p-value = 0.1085
## 
## Model df: 0.   Total lags used: 10
fc_ibmclose_nnetar <- forecast(ibmclose_nnetar, h=36)
autoplot(fc_ibmclose_nnetar)

The fitted value looks pretty good and residuals looks like white noise. However, even neural network method yielded naive-method like result in forecast. It looked like there wasn’t any rule in lagged values.

# experiment with usmelec data.
usmelec_nnetar <- nnetar(usmelec)
checkresiduals(usmelec_nnetar)

## 
##  Ljung-Box test
## 
## data:  Residuals from NNAR(3,1,2)[12]
## Q* = 215.96, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24
autoplot(usmelec) + autolayer(usmelec_nnetar$fitted)

fc_usmelec_nnetar <- forecast(
  usmelec_nnetar, h = 12*4
)

autoplot(fc_usmelec_nnetar)

It looked like neural network model was fitted well to the data.